
### This code summarizes MQ Supreme Court results
### Begins with posted results
### First, merges MCMCpack estimates against posted results, finds a correlation of r=0.99
##  Next, correlates VI estimates against MCMCpack results

rm(list = ls())
setwd("~/Documents/working/MCMC_mq")

load("MQout.Rdata")
mqdat <- read.csv("justices.csv")

#### Process MCMCpack output file (result)
result <- apply(out,2,mean)
lower <- apply(out,2,quantile, c(0.025))
upper <- apply(out,2,quantile, c(0.975))
lower <- lower[1:697]
upper <- upper[1:697]
result <- result[1:697]
tmp <- strsplit(names(result),"[.]t")
year <- length(tmp)
justice <- length(tmp)
for(i in 1:length(tmp)) year[i] <- as.numeric(tmp[[i]][2])
year = year + 1936
tmp <- strsplit(names(result),"[.]")
for(i in 1:length(tmp)) justice[i] <- tmp[[i]][2]
justice[which(justice=="O")] <- "O'Connor"

#### Merge main estimate from MCMCpack
mqdat$newest <- NA
for(i in 1:length(result)){
	index <- which(as.character(mqdat$justiceName) == justice[i] & mqdat$term == year[i])
	mqdat$newest[index] <- result[i]
}

### Correlation between MQ posted and MCMCpack-estimated results is 0.99
cor(mqdat$newest,mqdat$post_mn)

### Compare against VI results
load("mq_out.rda")

### Overall 0.93 correlation
vi.out <- c(t(lout$means$x))
vi.out[vi.out==0] <- NA
vi.out <- na.omit(vi.out)
cor(vi.out,result)	#0.93 overall
cor(vi.out[justice!="Douglas"], result[justice!="Douglas"])	#0.96 without Douglas

justiceNames <- unique(justice)
for(i in 1:length(justiceNames)){
  selector <- which(justice == justiceNames[i])
  cat("\n", justiceNames[i], ": ", signif(cor(vi.out[selector],result[selector]),digits=3), ", N=", length(selector),sep="")
}


### GRAPHICS: INDIVIDUAL LEGISLATORS OVER TIME

### Rescale VI results to MCMC scale
model <- lm(result ~ vi.out)
vi.out <- coef(model)[1] + coef(model)[2]*vi.out

## Merge bootstrap trials from VI estimates
load("bootstrap.rda")
vi.bootstrap <- matrix(NA, nrow=length(vi.out), ncol=length(vi.result))
for(i in 1:length(vi.result)){
	vi.trial <- c(t(vi.result[[i]]$means$x))
	vi.trial[vi.trial==0] <- NA
	vi.bootstrap[,i] <- na.omit(vi.trial)
}
## Rescale each iteration onto MCMC result 
for(i in 1:(length(vi.result))){
  model <- lm(result ~ vi.bootstrap[,i])
  vi.bootstrap[,i] <- coef(model)[1] + coef(model)[2]*vi.bootstrap[,i]
}
colnames(vi.bootstrap) <- paste("trial.", 1:ncol(vi.bootstrap), sep="")

## Individually Bias-correct the estimates
for(i in 1:nrow(vi.bootstrap)){
   bias = vi.out[i] - mean(vi.bootstrap[i,])
   vi.bootstrap[i,] = vi.bootstrap[i,] + bias
}

## Calculate the CIs
CI <- t(apply(vi.bootstrap,1,quantile, c(0.025,0.975)))
colnames(CI) <- c("lower","higher")
CI <- cbind(CI, mean=apply(vi.bootstrap,1,mean))


### PLOTS BY JUSTICE: FIGURE 9
## Excludes Byrnes, Sutherland, Cardozo, Brandeis, and Butler, who have less than 3 terms

for(i in c(1:32, 37:45)){

  pdf(paste("pdf/", justiceNames[i], ".pdf",sep=""))

  selector <- which(justice == justiceNames[i])
  title <- paste(justiceNames[i], " (r = ", signif(cor(vi.out[selector],result[selector]),digits=3), ", N=", length(selector), ")", sep="")
  plot(0,0,type="n",xlim=c(1930,2020),ylim=c(-6,6), xlab="", ylab="", main=title, bty="n", cex.main=2.3, cex.axis=2.0)
  polygon( c(year[selector],rev(year[selector])), c(lower[selector], rev(upper[selector])), col="grey", border=NA)
  lines(year[selector], vi.out[selector], lwd=5, lty=1, col="red")

##  If bias-correction by legislator is desired
#  bias = mean(vi.out[selector]) - mean(CI[selector,"mean"])
#  CI[selector,] = CI[selector,] + bias

  lines(year[selector], CI[selector,"lower"], lwd=2, lty=2, col="red")
  lines(year[selector], CI[selector,"higher"], lwd=2, lty=2, col="red")
  #lines(year[selector], result[selector], lwd=5, lty=3, col="black")

  if(i %in% which(justiceNames %in% c("Black","White"))) legend(1935, 6, lwd=c(5,4,15), lty=c(1,2,1), col=c("red","red", "grey"), legend=c("VI Estimate","VI 95% CI","Martin & Quinn 95% CI"), cex=2, bty="n")
  dev.off()
} #end for


sessions <- as.numeric(unlist(strsplit(names(result),"[.]t"))[seq(2,2*length(result),by=2)])

corr.r <- rep(NA, 77)
corr.spearman <- rep(NA, 77)
for(i in 1:77){
  corr.r[i] <- cor(result[sessions==i], vi.out[sessions==i])
  corr.spearman[i] <- cor(result[sessions==i], vi.out[sessions==i], method="spearman")
}

### PLOTS BY TERM: FIGURE 8
pdf("pdf/corr_by_term.pdf")
par(mar= c(5, 4, 10, 2) + 0.1)
plot(0,0,xlim=c(1930,2020), type="n", ylim=c(0.7,1), xlab="Year", ylab="Correlation")
points(1937:2013, corr.r, pch=1,  cex=1.2)
points(1937:2013, corr.spearman, pch=24 , bg="grey", cex=1.2)

legend(1985,0.9, pch=c(1,24), pt.bg=c("grey"), legend=c("Pearson","Rank Order (Spearman)"), bty="n")
dev.off()





